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ABSTRACT 

Two  simple  stationary  processes  of  discrete  random 
variables  with  arbitrarily  chosen  first-order  marginal 
distributions,  DARMA(p,N+l)  and  NDARMA(p,N),  are  given.   The 
correlation  structure  of  these  processes  mimics  that  of  the 
usual  linear  ARMA(p,q)  processes.   The  relationship  of  these 
processes  to  mover-stayer  models,  and  to  models  for  discrete 
time  series  given  separately  by  Lindqvist  and  Pegram  is  discussed 
Ad-hoc  nonparametric  estimators  for  the  parameters  in  the 
DARMA(p,N+l)  and  NDARMA(p,N)  are  given.   A  simulation  study 
shows  them  to  be  as  good  as  maximum  likelihood  estimators  for 
the  first-order  autoregressive  case,  and  to  be  much  simpler  to 
compute  than  the  maximum  likelihood  estimators. 
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1.   INTRODUCTION 
Discrete  time  series  arise  in  many  different  contexts. 
For  example  the  exact  arrival  times  in  an  arrival  process  are 
usually  not  measured.   Instead  the  number  of  arrivals  in  suc- 
cessive time  intervals  are  given.   This  is  the  case  with  the 
statistics  published  by  the  Center  for  Disease  Control  on  the 
incidence  of  various  diseases  in  the  United  States.   The  data 
are  given  as  the  number  of  occurrences  of  each  disease  in  suc- 
cessive days.   If  the  time  intervals  are  short  enough  and  the 
arrival  process  is  orderly,  then  the  resulting  time  series  is 
approximately  binary.   In  other  instances  the  process  that  is 
being  measured  is  continuous  but  the  data  is  quantized  in 
recording.   For  example,  the  amount  of  rainfall  in  a  day 
(24  hours)  at  a  location,  given  that  some  occurs,  is  a  continu- 
ous random  variable;  however,  it  is  often  recorded  to  the 
nearest  one-hundredth  or  one-tenth  of  an  inch.   Also,  since  a 
rainfall  series  will  often  contain  many  zeros  (no  rain) ,  an 
analysis  is  often  made  of  successive  wet  and  dry  days  which  is 
a  binary  time  series  [cf.  Buishand,  1978].   An  economic  impera- 
tive for  modelling  and  predicting  the  binary  rainfall  series  is 
that  it  is  the  primary  concomittant  variable  for  predicting 
volume  of  business  done  in  some  establishments  on  successive  days. 

Markov  chains  have  been  used  as  models  for  stationary 
discrete  time  series.   However,  they  are  overparametrized  for 
statistical  purposes.   Further,  the  data  to  be  modelled  can 
often  be  shown  to  be  non-Markovian,  or  at  least  not  first-order 
Markovian.   Higher  order  Markov  chains  can  be  used  but  this  only 
aggravates  the  problem  of  overparametrization. 


In  the  past  several  years  various  parametrically  simple 
models  have  been  proposed  for  stationary  discrete  time  series. 
The  models  have  as  parameters  the  fixed,  first-order  marginal 
distribution  of  the  time  series  and  the  correlation  structure. 
In  Jacobs  and  Lewis  [1978a,  1978b]  a  simple  scheme  is  given  for 
obtaining  a  stationary  sequence  of  discrete  random  variables  with 
a  given  marginal  probability  mass  function   tt   and  an  autocor- 
relation structure  like  that  of  a  mixed  first-order  autoregres- 
sive- (N+l) st-order  moving  average  process.   This  DARMA(1,  N+l) 
process  has  nonnegative  correlations  and  a  possibly  countably 
infinite  state  space.   The  correlation  structure  is  determined  by 
parameters  that  are  independent  of  the  marginal  distribution. 

A  special  case  of  the  DARMACl,  N+l)  process  with  mar- 
ginal probability  mass  function   tt   is  the  DAR(l)  process. 
This  is  a  Markov  chain  with  discrete  state  space  IE  and  with 
transition  matrix 

(1.1)  P  =  pi  +  d-p)Q  , 

where   Q   is  a  matrix  with   Q.  .  =  tt(j)   for   i,j£IE  ;   I   is  the 
identity  matrix  with  (i,j)  element   I. .   and   0  <  p  <  1  .   The 
correlation  structure  of  a  real  valued  DAR(l)  i.e.  one  for  which 
IE  is  a  subset  of  the  real  line,  is  that  of  a  first-order  auto- 
regressive  process  with  kth-order  serial  correlation  equal  to 
p   .   There  is  no  limitation  on   tt  ;  a  common  and  useful  assump- 
tion is  that  it  be  Poisson  and  therefore  have  an  infinite  state 
space.   The  DAR(l)  model  with  a  finite  state  space  is  a  special 
case  of  the  mover-stayer  model  [Bhat,  1972,  p.  302-9]. 

Lindqvist  [1978]  studied  a  real  valued  finite  state 
space  Markov  chain  with  a  transition  function  of  the  form 


(1.2)  P, 


ij  =  CIij  +  d-d  Qi 


where   Q.  .  =  tt  ( j )  ,  for   i,  j  6  IE  as  before.   Since  the  state  space 
IE  is  finite,  the  constant   c   can  take  on  some  negative 
values  with  the  constraint  that 


(1.3)  max  [1  -  (l-Ti(i))  X]  <  c  <  1  , 

l<i<r 


where   r   is  the  number  of  elements  in  the  state  space. 

In  Jacobs  and  Lewis  [1978c],  the  DAR(l)  process  was 
extended  to  obtain  a  sequence  of  discrete  random  variables 
with  pth  order  Markov  dependence  and  given  marginal  distri- 
bution.  The  DAR(p)  process  is  defined  as  follows.   Let  (V  } 
be  a  sequence  of  independent  identically  distributed  random 
variables  with   P{V   =  1}  =  1  -  P{V   =  0}  =  1  -  p, 
0  <  p  <  1;  {A  }   is  a  sequence  of  independent  identically 
distributed  random  variables  taking  values  {l,2,...p}  with 
P{A   =  i}  =  a./  i  =  1 , 2 , . . . , p ;  and  {Y  }  is  a  sequence  of 
independent  identically  distributed  random  variables  with 
discrete  state  space   IE  and   P{Y   =  i}  =  tt  ( i )  .   Let 


(1.4)  Z   =  V  Z   A   +  (1-V  )Y 

n     n  n-A         n   n 
n 


The  process  {Z  }  is  called  a  DAR(p)  process.   Note  that  by 
direct  argument  from  (1.4) 


(1.5)       P{Z  ...  =  jlz   =  i,,...,Z    ,,  =  i  } 
n+1    Jl  n     1'    '  n-p+1     p 


P 
(1-p)  tt(j)  +  I      pa  e  (i  )  , 
k=l    K  J      k 


where   e  .  (i)  =  1   if   i  =  j   and   e  .  (i)  =  0   otherwise.   There 
is  no  limitation  on  the  marginal  probability  mass  function   tt  . 

Pegram  [1980]  considers  a  real-valued  finite  state  space 
model  {Z1}  which  is  a  generalization  of  the  DAR(p)  model  in  that 
its  conditional  probabilities  are  of  the  form 


d'6)       P{Zn+l  =  ^Zn  =  h Zn-P+1  =  V 


k=l     k  k=l     K  J  K 


where   Hi,;  k=l,...,p}   are  (possibly)  negative  constants.   Note 

that  although  some  of  the  constants   (}>,   may  be  negative,  the 

admissible  values  for   {i  }   depends  on  the  marginal  distribution 

tt  .   It  was  shown  in  Jacobs  and  Lewis,  [1978c]  that 

Corr(Z  ,  Z  ,,),  k=l,2,...   for  the  real  valued  DAR(p)  process 
n    n+k  ^   ^ 

are  nonnegative.   Pegram' s  model  allows  some  of  the  correlations 
to  be  negative.   The  amount  of  negative  correlation,  as  in 
Lindqvist's  model,  depends  on  the  marginal  distribution   tt  . 

In  this  paper  we  will  consider  models  for  real-valued 
stationary  discrete  time  series  whose  nonnegative  correlation 
structure  is  that  of  a  mixed   pth-order  autoregressive  and   qth-order 
moving  average  process.   Thus  we  have  a  generalization  of  all 


of  the  preceding  models.   In  Section  2  we  will  give  definitions 
of  two  such  models,  DARMA(p,N+l)  (discrete  mixed  autoregressive- 
moving  average  process  with  orders   p   and   N+l  respectively) 
and  NDARMA(p,N).  We  briefly  describe  some  of  their  properties 
and  suggest  an  estimate  for  the  correlations.   In  Section  3 
we  describe  in  detail  a  simulation  experiment  that  was  done 
to  study  the  behavior  of  various  estimators  for  the  first  order 
serial  correlation  coefficients   p   of  the  DAR(l)  model  for 
small  and  moderate  sample  sizes.   In  Section  4  some  extensions 
of  the  DARMA  models  are  briefly  discussed  including  one  which 
can  have  negative  correlations.   Throughout  the  remainder  of 
the  paper  we  will  assume  that  the  NDARMA  and  DARMA  processes 
are  real  valued.   They  can  in  fact  be  used  to  model  categorical 
time  series,  but  then  numerical  measures  such  as  correlations 
are  meaningless. 


2.   The  DARMA(p,N+l)  and  NDARMA(p,N)  Processes. 
In  what  follows  we  let  {Y  }  be  a  sequence  of  independ- 
ent identically  distributed  random  variables  taking  values  in 
a  real-valued  discrete  state  space   IE   with  P{Y   =  i}  =  tt  ( i)  , 
ielE .   Let  {U  }  and  {V  }  be  independent  sequences  of  independ- 
ent random  variables  taking  the  values   0   and   1   with 


(2.1)       P{Un  =  1}  =  e     and     P{Vr  =  1} 


for  fixed   0    $    1   and   0  £  p  <  1.   Let  {D  }  be  a  sequence 

of  independent  identically  distributed  random  variables  taking 

values  0,1,2, .  .  .  ,N  with  P{D,  =  n}  =  6   ,  n  =  0,1,. ..,N,  and 

x.  n 

{A  }  be  a  sequence  of  independent  identically  distributed 
random  variables  taking  values  l,2,...,p  with  P{AV    n}    an  , 
n  =  1,2, ... ,p. 


2.1   The  DARMA(p,N+l)  process. 

The  DARMA(p,N+l)  process  is  a  sequence  of  random 
variables  {X  }  which  is  formed  according  to  the  probabilistic 
linear  model 


(2.2)         X   =UY   _   +  (1  -  U  )Z   ,„  -,  . 
n     n  n-D  n   n-(N+l) 


for   n  =  1,2,...,  where  the  "autoregressive  tail"  is 


(2.3)  Z   =VZ   ^   +  (1  -  V  )Y 

n     n  n-A  n   n 

n 


for   n  =  -N-p+1,  -N-p+2,....   This  process  differs 
from  the  DARMA(1,N+1)  process  defined  in  Jacobs  and  Lewis 
[1978a]  in  that  the  "autoregressive  tail"   Z    is  now  the   pth 
order  autoregressive  process  of  (1.4).   In  Jacobs  and  Lewis 
[1978c],  it  was  shown  that  the  vector-valued  Markov  chain 
{(Zn,  zn+1'---'  zn_p+1)  '  n  =  1,2,...}  has  a  limiting  joint 
probability  mass  function   v   with  marginal  probability  mass 
function  i\ .   Hence,  if  (Z_N_   ,  ,  .  .  .  ,  Z_  )  has  joint  probability 
mass  function   v  ,  then  {X  ;  n  =  1,2,...}  is  a  stationary 
process  with  marginal  probability  mass  function   tt  . 

Let   r(k)  =  Corr(X  ,  X  ,,  )  for  the  stationary  process. 

n    n+k  J    c 

Then  {r(k)}  can  be  shown  to  satisfy  the  following  system  of 
equations : 

7    N_1  ? 

(2.4)     r(l)  =  3    I      <5.6     +  3(l-3)rR(l)  +  (1-3)  r.(l)  , 

i=0   x  x  x  a  A 


N-2 
(2.5)     r(2)  =  3    I   6i6i  +  2  +  B(l-3)rB(2)  +  (1-3)  r  (2)  , 

i  =  0 


(2.6)     r(N)  =  326Q6N  +  3(l-3)rB(N)  +  (l-3)2rA(N)  , 


(2.7)     r(N+k)  =  3(1-3)  r_(N+k)  +  ( 1-3)  2rA  (N+k)  .      k    1 


In  these  equations 


(2.8)     rA(k)  =  Corr(Zn,  ZR+k)  ,  k  >  1, 


which  satisfy  the  following  Yule-Walker  equations: 


(2.9)     rA{1)    =    PairA(°)  +  Pa2rA(1)  +  •••  +  paprA(p_1)  ' 


(2.10)  rA(2)  =  pairA(l)  +  pa2rA(0)  +  ...  +  pa  rA(p-2)  , 

•  •  •  • 

•  •  •  • 

(2.11)  rA<P)  =  Pa1rA(P-D  +  pa2rA(P~2)  +  ...  +  pa  rA(0)  , 

and  for   k  :>  1  , 

(2.12)  rA(p+k)  =  pa1rA(p+k-l)  +  pa2rA(p+k-2)  +  . . .  +  pa  r  (k)  , 

where   r  (0)  =  1 . 
In  addition 


rB(i)  =  Corr(Zn+ 


N 


i-(N+l)'  Yn-Dn)  =  lQ    Corr(Zn+i-(N+l)'  Yn-j)6j 


is  obtained  recursively  as 


rB(0)  =  0  ; 


rB(l)  =  (1-P)6N  ; 

rB(2)  =  POlrB(l)  +  d-P)Vl 


k-1 
rB(k)  =  J^  pa.rB(k-i)  +  (1-p)  6N_(k_1} 

for   k  <  min(p,N)  ; 


rfi(k)  =  Ji  pa.rB(k-i)  +  (1-p)  6N_(k_1} 


for  max(p,N)  >  k  >_  min(p,N)  • 


rfi(k)  =  l      pairB(k-i) 
i=l 


for   k  >  max(p,N)  . 


To  see  that  the  serial  correlations  for  the  DARMA (p, n+1) 
process  are  all  nonnegative  let   q(i)   (respectively  q  (i))  be 

the  probability  that   X    and   X      (respectively   Z    and 

^  n        n+i      ^  J         n 

Z  ,.)   choose  the  same  random  variable   Y.  ,  where,  because  of 
n+i  k 

the  backward  definition  of  the  autoregression  k  <  n.   Then 
q(i)  (respectively  q  (i))  also  satisfy  equations  (2.4)  -  (2.7) 
(respectively  (2.9)  -  (2.12))  and  since  they  are  nonnegative, 
so  are  the  serial  correlations. 

To  see  this  identity,  let   R   be  the  random  index  of 

1  n 

the   Y,  ,  k  <  n,   that   X    chooses;  that  is, 
k '    —  n 


X   =  YD 
n     R 
n 


Then,  since  the  random  variables   R    are  independent  of  the 

n 

fr,  }  random  variables, 


(2.13)     E 


[x  x  _,.]  =  e[y_  yd 

n  n+£        R    R  ,0 
L       J      |_  n    n+£ 


n 

I       E 
k=l 

S 

P 

n      n+£ 

+   l     l 

E 

k=l   j=l 

J7* 

P{R   =  R     =  k} 
n     n+x, 


P{R   =  k,  R  ^0  =  j} 
n        n+£    J 


(2.14) 


=  E 


K 


P{R   =  R  .  . , 
n     n+£j 


Thus 


+  E 


P{Rn  *  Rn+£} 


Cov(X  ,  X    ) 

n    n  +  x, 


=  E 


P{R   =  R    } 
n     n  +  £ 


+  E 


3[Yl]   [P{Rn  *   Rn  +  £>  -  l 


(2.15) 


H 


Vai  |Y, I  P{R   -  R  .  „} 

n     n+£ 


Therefore 

(2.16)     Corr(Xn,  Xn+S,)  =  P{Rn  =  Rn+£}  -  q(£) 

as  asserted  above. 

This  identity  will  also  be  used  in  the  estimation 
procedure  proposed  in  Section  3  for  the  serial  correlations 


10 


2.2   The  NDARMA(p,N)  process 

In  this  subsection  we  will  define  another  related 
discrete  time  series  with  the  correlation  structure  of  a  mixed 
moving  average  autoregressive  process.   This  new  process  is 
more  reminiscent  of  the  linear  ARMA(p,N)  process.   The  key  idea 
that  leads  to  this  new  model  is  that  a  probabilistic  mixture  of 
a  finite  number  of  random  variables  each  with  probability  mass 
function   tt   has  probability  mass  function   it   even  if  the 
random  variables  are  dependent.   Thus  it  is  not  necessary  to 
define  the  autoregression  via  an  autoregressive  tail,  as  in  the 
DARMA(p,N+l)  process;  the  autoregression  can  be  made  explicit, 
as  in  the  usual  (normal  theory)  linear  processes. 

Thus  let 


<2-17'  Xn  "  VnXn-A   +  (1-Vn)Yn-Dn 

n  n 


where  {V  } ,  {A  } ,  and  {D  }  are  as  before.   Thus,  with  prob- 
ability  p  ,   X    is  one  of  the   p   previous  values 

X   ,,..., X     and  with  probability  (1-p)  it  is  a  mixture  of 
n-1       n-p 

the  previous   Y's,n-N<k£n.   Note  that  if   p  =  0, 

then   {X  ;  n  =  1,2,  .  . .  }   is  a  DMA(N)  process  as  defined  in 
n 

Jacobs  and  Lewis  [1978a].   If   P{Dn  =  0}  =  1,  then   {Xr}   is 

a   DAR(p)  process  as  defined  in  Jacobs  and  Lewis  [1978c]. 

Let   t  =  inf{i  :  6.  >  0}  .   Note  that 

Z   =i(X    ,X   t  ,   ....  X      -.   ,   x      ,...,Y   ,.  /   ,  n  •   X  ,_,...  i 
^n       n  '   n-1         n-p+1     n-T         n-N 

is  a  Markov  Chain  with  state  space  IF   which  is  equal  to  the 
product  space  of   IE  with  itself   p  +  (N-x)   times.   Since 
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P{Xn+T+l   =    Yn+1   =    ^X0 Xn'Y0 Yn}    >    (1"p)(St    *  (  V  > 


there    is    a    set      JcrlF       such    that      min   P(Z      „    =   kjz      =   4J    =   y    >    0    , 

4^IF    n+K      ^ 

where 

K  =  P  +  N   and 

{Xn+K    Yn+K-x'  Xn+K-1    Yn+K-l-x ' * *  * '  Xn+K-p+l  =  Yn+K-p+l-x} 

Thus  the  condition  of  case  (b)  on  page  173  of  Doob  [1953]  is 
satisfied.   The  proof  on  pagesl73  and  174  extended  to  countable 
state  spaces  shows  that   &   has  a  limiting  probability  mass 
function   v   as   n  ■*  °°  ;  further  the  convergence  of  the  con- 
ditional distribution  of   Z    to   v   as   n  -*■  °°   is  geometric. 

^n  3 

The  marginal  probability  mass  function  of   v   is   tt  . 

It  follows  from  (2.17)  that  the  serial  correlations 
for  the  stationary  NDARMA(p,N)  process  satisfy  the  Yule-Walker 
equations  for  the  ARMA(p,N)  process  with  restrictions  on  the 
range  of  the  coefficients; 


(2.19)       rN(k)  =  Corr(Xn,  Xn+k) 


=  .^  P\    Corr(Xn,  Xn+k_.) 

+  (1-p)  J^  «iCorr(Xn  ,  Yn+k_.) 
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for   k  >  0  .   The  correlations   rMD(i)  =  Corr(X   ,  Y   .)   can 
—  NB  n     n-i 

be  computed  recursively  as  follows.   For   i  =  0 


(2.20)  rNB(0)  =  Cl-P)6o  ; 


for   1  <  i  <  p 


rNB(i)  =  (1"P)6i  +  PalrNB(i~1)  +  •••  +  Pairnb(0)  '" 


and  for   i  >  p 


rNBU)  =  d-P^i  +   X   PairNBU-^   ' 


where  if   i  >  N  ,  then   6.  =  0   by  convention. 
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Hence,  if  we  assume   N  <  p 


(2.21)     r  (k)  =  pa1rN(k-l)  +  pa2rN(k-2)  +  ...  +  paprN(p-k) 


+  Cl-P)   f   rNB(i"k) 
i=k 


for   1  <  k  <  N  ; 

(2.22)     rN(k)  =  pa1rN(k-l)  +  pa2rN(k-2)  +  ...  +  pa  rN(k-p) 

for   k  >  N  . 

The  serial  correlations  of  the  NDARMA(p,N)  process  are 
nonnegative  since,  if   ^(i)   i-s  the  probability  that   X    and 

X  , •   choose  the  same  random  variable   Y,  ,   k  <  n  ,  then 

n+i  k       —    ' 

{q  (i)  }  satisfies  equations  (2.21)  and  (2.22).   The  argument 
is  the  same  as  for  the  DARMA(p,N+l)  case. 

2.3   Comparison  of  Admissible  Range  of  Correlations  for  the 
DARMA(1,1)  and  the  NDARMA(1,1). 

Let  (Xn)  be  a  stationary  DARMA(1,1)  process;  that  is, 


(2.23) 


X   =  < 
n 


n 


'n-1 


with  probability  3, 
with  probability  1-3, 
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where 
(2.24) 


Z   =  i 
n 


'n-1 


n 


with  probability  p, 
with  probability  1-p. 


Equations  (2.4)  -  (2.12)  for  the  DARMA(p,N+l)  correlations 
simplify  to 


(2.25)     r(k)  =  Corr(Xn,  Xn+R)  =  pk  1  (1-3)  [ B  (1-p)  +  (l-B)p] 


Similarly  let  {x'}  be  a  stationary  NDARMA(1,1)  process; 
that  is 


(2.26) 


X* 
n 


n-1 


=  < 


n 


Y 


n-1 


with  probability  p, 


with  probability  (l-p)6n, 
with  probability  (1-p) (1-6-J 


Equations  (2.21)  and  (2.22)  simplify  to 


(2.27) 


rM(k)  =  Corr (X' ,  X '  ,  ) 
N  n    n+k 


k-1 


[p  +  (l-p)260(l-60)  1 


Figure  1  gives  graphs  of  the  attainable  values  of 
{r  (2) ,  r  (1) }  as  the  parameter  values   p  and  6_   vary,  and 
{r(2),  r(l)}  as  the  parameters   p  and  3  vary.   Note  that 
although  the  set  of  attainable  correlations  for  the  NDARMA(1,1) 
process  is  not  strictly  contained  in  that  for  the  DARMA(1,1) 
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process,  it  is  much  smaller.   Thus,  the  DARMA(1,1)  model  appears 
to  be  broader  than  the  NDARMA(1,1)  model.   The  smaller  region 
of  possible  correlation  pairs  for  the  NDARMA(1,1)  model  seems  to 
be  a  constraint  due  to  the  explicit  autoregression  on   X  _,  . 

2.4   An  Estimator  for  the  Serial  Correlations  of  the 
DARMA( p,N+l)  and  NDARMA(p,N)  processes. 

The  usual  estimator  for  the  serial  correlations  of  a 
stationary  real-valued  sequence  {X,  ,  .  .  .  ,  X  }  is 


(2.28)      r(Z)     =    [S*]  x  (m-i)     x  J       (X.  -  X) (X.+l    -  X)  , 


where 

-1   m 
(2.29)  X  =  m     I      X. 

j=l  3 


and 


2     1    ?   ,„     ;,2 


(2-30)  ^-   i^T  .1       <Xj 


In  this  subsection  we  will  suggest  another  estimator 
for  the  serial  correlations  of  the  DARMA  and  NDARMA  processes 

By  the  remarks  at  the  ends  of  Sections  2.1  and  2.2, 
the  £th  serial  correlation,  r(£),  for  both  the  stationary 
DARMA(p,N+l)  and  NDARMA(p,N)  processes  is  equal  to  the  prob- 
ability that   X    and  *n+i      choose  the  same   YR  ,  k  <  n  . 
Hence,  for  both  processes  for   i  f    j, 
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P{Xn  =  i,  X     =  j}  =  P{Y    =  i,  Y 

n        n+£  Rn        Rn+, 


=  j) 


I  P(Rn  =  k,  Rn+£  =  rl  P{Yk  =  i,  Yr  =  j 

k<n  r£n  +  £ 


Since  the  {Y,  }  random  variables  are  independent,  and  independent 

K. 

of  R   and  R  , 0  ,  we  have 
n       n+£ 


(2.31)       P{Xn  =  i,   Xn  +  £  =  j}  =  TT(i)TT(j)P{Rn  ft     Rn  +  £} 

=  7T(i)7T(j)  [1  -  r(£)  ]   , 

by  equation  (2.16). 

Thus,  for   j  €  IE 

,  N-m 

(2.32)       lim   B    <m,j)  lim    (N-m)  £         £       1     (x    )     1    (X,        ) 

N^oo  N->oo  i^j    k=i      !      K         :      k+m 

=    [1-tt(j)  ][l-r(m)  ]    7T(j) 


almost  surely  where   1± (x)  =  1   if   x  =  i   and   0   otherwise. 
Hence 


(2.33)       r(m)  =  1  -  I         B  (m,j)[l  -  tt  (  j  )  ]  1 

j  €IE 


is  a  strongly  consistent  estimator  for  the  mth  serial  correla- 
tion for  the  stationary  DARMA  and  NDARMA  processes. 


Estimator   r (m)   is  also  a  strongly  consistent  estimator  for 
the  finite  state  space  models  of  Lindqvist (1978)  and  Pegram 
(1980)  since  the  conditional  probabilities  (1.2)  and  (1.6) 
are  of  the  same  form  as  those  for  the  appropriate  DAR(p) 
process . 

In  the  next  section  we  pursue  this  estimator  for  the 
special  case  of  the  first-order  autoregressive  process  DAR(li 
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3.   ESTIMATION  FOR  THE  DAR(l)  PROCESS 
Let   {X  }   be  a  stationary  DAR(l)  process  with  state 
space   IE  =  {0,1,...}   and  first-order  serial  correlation   p  , 
0  <  p  <  1  ;  that  is, 


(3.1)      Xn  =  i 


X  _,  with  probability  p  , 


Y  with  probability  1-p 


for   n  =  1,2,...,  while   Xn   is  a  random  variable  independent 
of   {Yn}   but  with  the  same  probability  mass  function. 

A  simulation  was  conducted  to  study  the  performance  of 
several  estimator  for   p   for  small  to  moderate  series  lengths 
m  .   The  series  lengths  considered  were   m  =  20 ,  50,  and  200. 
The  marginal  probability  mass  functions  considered  were  the 
Poisson  with  parameter   A  , 


-A  Ak 
(3.2)        TT(k)  -  e   £y        k  =  0,1, 


and  the  geometric  with  parameter   p 


(3.3)         7T(k)  =  pk(l-p)        k  =  0,1, 


One  type  of  estimator  considered  was  the  single  param- 
eter maximum  likelihood  estimator.   For  a  series  of  length   m 
let   N. .   denote  the  number  of  times  the  DAR(l)  process  goes 
from   i   to   j  ,  for   i,  jelE   and  let   N.   denote  the  total 
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number  of  times  the  process  is  in  state   j  .   The  log-likelihood 

function  for  a  DAR(l)  series   {X,,...,X  }   of  lenqth   m   is 

l       m  3 


(3.4) 


=  I       J      Ni   £n[(l-p)Tr(j)  ] 


i=0  jfti 


+      I      N.  .  £n[l  -  (l-p)TT(i)  ] 
i  =  0 


+   7   1.  (X,)  tt  (i)  . 
i=0  x      x 


Taking  the  partial  derivative  of   L   with  respect  to   x  =  1-p 
and  setting  the  derivative  equal  to  zero  results,  after  some 
simplification,  in  the  following  equation  for  the  maximum  likeli- 
hood  estimator   x  =  1-p  ,  if  it  exists: 


(3.5)     f(x)  =    1  -  N     I      N..  {1  -  x[l  -  ir(i)]}    =  0  . 

i=0   1:L 


Note  that   f(x)   is  monotone  decreasing  in   x   and   f(0)  >_   0 
Hence,  if  there  is  a  solution  to  (3.5)  in  [0,1],  it  will  be 
unique. 

The  ad-hoc  estimator  of   p   given  at  (2.33)  was  also 
considered.   For  the  DAR(l)  process  this  estimator  is 


pi  [N  1       I       N.   ][1  -  7T(j)  ]  1 

a  — n  i       ■;  — n   1-J 


1 


j=0 


i  =  0 


I    * 


21 


The  summations  in  both  the  maximum  likelihood  estimators 
and  the  estimators  of  the  form  (3.6)  for  the  Poisson  case 
(respectively  the  geometric  case)  were  restricted  to  be  between 
a  =  max[l,  n_(y-7o)]   (respectively   a  =  max ( 1 ,  n_(y-10a)) 
and   b  =  n  (y+7a)  (respectively  n+ ( y+10o) ) ; here  n_(y)  (respective 
n  (y) )  is  the  largest  (respectively  smallest)  integer  less 
(respectively  greater)  than   y  . 

Equation  (3.5)  was  solved  numerically.   In  the  case 
N  =  20  ,  it  was  not  uncommon  that   f (x)   did  not  have  a  zero  in 
[0,1].   In  this  case,  if   f(l)  >  0  ,  then   x   was  set  equal  to 
1  ;  that  is,   p=0.   If   f(0)=0,  then   x   was  taken  to  be 
0  ;  that  is,   p  =  1  . 

Other  estimators  for   p   that  were  considered  included 
the  following: 

1.   The  usual  estimate  for  first-order  serial  correlation 


(3.7)     p±    =    [s£]    (m-1)  L       I       (X   -  X) (X     -  X) 

n=l 


-         2  2 

where   X   and   S    are  as  in  (2.29)  and  (2.30) .   If   Sx  =  0  , 

then   p,   was  set  equal  to   1  . 

2.   The  maximum  likelihood  equation  (3.5)  was  solved  numer- 
ically for   p   for  each  of  the  following  three  values  for   tt  ( j ) 

a.  the  known  distribution  (3.2)  or  (3.3)  with  known 
parameter   A   or   p   was  used;   p2   denotes  this  estimator; 

b.  the  known  distribution  with  an  estimated  parameter 
was  used;  that  is, 
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"■  k 

(3.8)  7T(k)  =  e"A  ^r-  where       A  =  X 

or 

(3.9)  TT(k)  =  pk  (1-p)       where       p  =  X[l  +  X]"1 


P-.   denotes  this  estimator; 

-1 

c.  the  nonparametric  estimator   tt(j)  =  N    N.   of 

tt(j)   was  used;   p.   denotes  this  estimator. 

d.  p^   is  the  estimator  of   P   resulting  from  the  two- 
dimensional  maximum  likelihood  estimator  where  the  other 
parameter  is  the  distribution  parameter  (A  or  p) ; 

3.  the  estimate   pfi   is  the  nonparametric  estimate 
(3.6)  using   N    N  .   as  the  estimate  for   tt  ( j )  ; 

4.  the  estimate   p7   is  the  nonparametric  estimate 
(3.6)  using  the  true  value  of   tt  ( j )  • 

Both  estimators   pfi   and   p7   can  have  negative  values 
for  small  to  moderate  sample  sizes.   Hence,  we  also  considered  the 
following  estimate. 

5.  Estimator   p0  =  max(p,  ,  0  ) . 

o  b 


3.2   The  sampling  experiment. 

A  DAR(l)  series  of  length   m  was  simulated  and  the 
estimates  for   p   were  computed.   The  computation  was  repeated 
for  1000  independent  replications  and  the  sample  mean,  sample 
variance,  and  sample  root  mean  square  error  were  computed.   Each 
experiment  was  then  repeated  for  20  independent  replications, 
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and  the  mean  of  the  means,  mean  of  the  standard  deviations,  and 
mean  of  the  root  mean  square  errors  were  computed.   Tables  1-3 
give  the  means  of  the  root  mean  square  errors  for  the  cases 
studied.   The  box  plots  of  row  values  appearing  in  the  last 
column  of  the  table  are  given  to  help  the  reader  to  summarize 
the  performance  of  the  8  estimators  across  the  7  cases  considered 

All  runs  were  performed  on  an  IBM  system  360/67  computer 
at  the  Naval  Postgraduate  School  using  the  LLRANDOM  package 
[Learmonth  and  Lewis,  19  73]  which  generates  numbers  according 
to  the  scheme  given  by  Lewis,  Goodman,  and  Miller  [1969].   Tests 
of  the  random  number  generator  are  given  in  Learmonth  and  Lewis 
[1974]. 

Among  all  the  estimates  of   p  ,  the  usual  first-order 
serial  correlation  estimator,   p,  ,  performs  least  well.   The 
maximum  likelihood  estimator  with  smallest  root  mean  square 
error  tends  to  be   p~  ,  although  by  the  time   m  =  200   the 
difference  is  minor.   Of  course  the  value  of   p   or   A   would 
not  be  known  in  general,  so  that  this  estimator  is  unrealistic. 
Maximum  likelihood  estimators   p_   and   p^   are  about  equiva- 
lent, indicating  that  the  extra  computational  complexity  of 
the  two-parameter  maximum  likelihood  estimator,   p^  ,  is  not 
necessary.   The  performance  of  the  nonparametric  estimator   p0 

o 
is  about  the  same  or  sometimes  better  than  that  of  the  maximum 

likelihood  estimator   p4  ,  especially  if,  for  small  sample  size 
(m=20) ,  the  modification   pg  =  max ( p  ,0)   is  used.  This  sug- 
gests that  the  ad-hoc  estimator  of  (3.6)  altered  to  give  values 
in  the  range  [0,1]  is  almost  as  good  as  the  maximum  likelihood 

24 


estimator  with  the  same  estimate  for   tt  ( j )  .   The  ad-hoc 
estimator  is  much  easier  to  compute  than  the  maximum  likelihood 
estimator. 
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4.   EXTENSIONS 

Since  the  process   {X  }   is  obtained  as  a  probabilistic 

mixture  of  the   {Y  },  the  DARMA  and  NDARMA  process  may  be  defined 

using  any  sequence  of  independent  identically  distributed  random 

variables   {Y  }  .   One  implication  is  that  DARMA  and  NDARMA 
n 

processes  may  have  a  continuous  marginal  distribution.   However, 
even  if  the  distribution  of   Y    is  continuous,  a  realization 
of  the  sequence   {X  }   will,  in  general,  contain  many  runs  of  a 
single  value.   This  seems  to  be  the  major  drawback  to  using 
DARMA  and  NDARMA  processes  to  obtain  a  sequence  of  dependent 
random  variables  with  a  specified  continuous  distribution  and 
correlation  structure.   However,  the  process  with  continuous 
marginals  may  be  useful  in  simulation  studies. 

Multivariate  DARMA  and  NDARMA  processes  may  be  obtained 
by  using  a  sequence  of  multivariate   Y  's  .   To  illustrate  this 
we  generate  DARMA  and  NDARMA-like  processes  having  negative 
correlations.   These  can  be  derived  from  bivariate  processes 
as  follows. 

Let   { (Y  (1),  Y  (-1))}   be  a  sequence  of  independent 
bivariate  random  variables  with  state  space   IE  =  [0 ,+!,... J     , 
marginal  probability  mass  function   tt  ,  and  correlation 
r  =  Corr(Y  (1),  Y  (-1))   which  will  be  negative  in  general.   One 
way  to  generate  such  a  sequence  is  to  note  that  a  random  vari- 
able  Y  (1)   with  probability  mass  function   tt   can  be  simulated 
from  a  uniform   [0,1]   random  variable  by  defining 


j-1  J 

(4.1)     Y  (1)  =  j       if        I       tt  (i)  <  U  <  I       TT(i) 
n  i=_oo  i  =  -<» 
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If   Y  (-1)   is  generated  by 


j-1  3 

(4.2)      Y  (-1)  =  j        if         I       TT(i)  <  1  -  U  <  I       7T(i)   , 

n  j_  =  _co  j_=— °° 


then   (Y  (1),  Y  (-1))   is  called  an  antithetic  pair.   If   tt   is 
n       n 

symmetric  about  zero,  then   Y  (1)  =  -  Y  (-1)  ,  and 
r  =  Corr(Yn(l),  Yn(-1))  =  -1  . 

A  bivariate  DARMA  (p,N+l)  process   {(X  (1),  X  (-1))} 
is  defined  as  follows.   Let   {an,...,a  }   and   {bn,...,b.-}   be 
fixed  sequences  of  numbers  that  are  either   -1   or   1  .   Let 


(4.3)      Xn(l)  -UnY     Cb   >  +  (1-VZn-(N+l)   (V 


(4.4)     Xn(-1)  =  UnYn_D  (-bD  )  +  d-On)Zn_(N+    (-a0) 

n     n 


for   n  =  1,2, 
where 


(4.5)     Zn(a0)  =  VnZn_A  (aA  )  +  d-Vn)Yn(a0 

n    n 


(4.6)     Zn(-a0)  =  Vnz     (-aA  )  +  (l-vn)Yn(-a0) 

n     n 


for   n  =  -N,  -N+l,...  where   {A  }   and   {D  }   are  as  in  Section 

n  n 

2.   The  random  variable   X  (-1)   is  called  the  dual  of   X  (1) 

n  n 
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or  the  antithetic  when  (4.1)  and  (4.2)  hold  for  the   {Y  }   pair 


n 


Note  that  if  (4.1)   and   (4.2)   hold  and   it   is  symmetric  about 

zero,  then   Z  (-1)  =  -  Z  (1)   and   X  (-1)  =  -  X  (1)  . 
n  n  n  n 

•       t 
A  bivariate  NDARMA(p,q)  process   (X  (1),  X  (-1))   can 

be  defined  similarly: 


(4.7)     Xn(l)  =  VnXA.fl  (aA  )  +tt-Vn)Yn_D  (bD  )  , 

n    n  n    n 


(4.8)    X;(-l)  =  vnx;_ft  (-a  )  +  d-vn)Yn_D  (-b   ) 

n     n  n     n 


The  stationary  bivariate  DARMA  and  NDARMA  processes 
will  have  marginal  probability  mass  function   tt  .   A  process 
having  possibly  negative  correlations  can  be  obtained  by  consider 
the  marginal  processes   (X  (1)},  {X  (-1)},  {X  (1)},  {X  (-1)1  . 
Details  will  be  given  elsewhere. 
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